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ABSTRACT 


Stenoses  (blockages)  in  coronary  arteries  cause  heart  attack,  which  the  Centers  for 
Disease  Control  and  Prevention  ranks  as  the  leading  cause  of  death  among  American 
adults.  We  have  performed  computer  simulations  in  order  to  understand  the  flow  patterns 
due  to  stenoses  with  and  without  a  guidewire  used  for  interventional  procedures  (e.g., 
stent  deployment).  Building  off  previous  models  that  have  been  partially  validated  with 
experimental  data,  this  thesis  continues  to  develop  the  models  in  order  to  further  evaluate 
fluid  characteristics  for  cardiovascular  applications.  Reasonable  agreement  occurs 
between  computational  and  experimental  data.  This  work  is  done  in  preparation  for 
future  research  to  develop  a  microelectromechanical  system  capable  of  monitoring  blood 
flow  and  aiding  in  medical  diagnoses. 
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I.  INTRODUCTION 


A.  BACKGROUND 

The  leading  cause  of  death  in  United  States  adults  is  heart  attack,  according  to  the 
Centers  for  Disease  Control  and  Prevention.  This  is  generally  caused  by  a  constriction  of 
a  coronary  artery  through  a  process  known  as  atherosclerosis,  or  an  increased  amount  of 
plaque  building  on  the  walls  of  the  artery  eventually  restricting  flow.  This  can  be 
particularly  dangerous  due  to  the  risk  of  a  blood  clot  (known  as  thrombus)  forming  on  the 
surface  of  the  plaque  or  a  piece  of  the  plaque  breaking  off.  Either  of  these  situations 
could  result  in  a  complete  blockage  of  an  artery  resulting  in  life  threating  medical 
conditions  (e.g.,  heart  attack). 

The  most  common  invasive  treatment  for  a  blocked  artery  is  the  use  of  a 
technique  known  as  angioplasty.  In  this  procedure,  a  balloon  is  carefully  placed  within 
the  artery  and  inflated  in  order  to  re-widen  the  artery  and  increase  blood  flow.  Generally 
a  stent,  or  a  mesh  wire  tube,  is  permanently  inserted  in  order  to  scaffold  open  the  artery 
wall  following  the  angioplasty  (Denardo  et  ah,  1994). 

Unfortunately,  stents  have  a  non-trivial  rate  of  failure.  Depending  on  patient 
characteristics  (e.g.,  diabetes),  artery  characteristics  (e.g.,  diameter)  and  design  of  the 
stent  (e.g.,  bare  metal  vs.  drug  eluting),  the  risk  of  the  stenosis  (blockage)  returning  and 
causing  failure  can  be  as  high  as  30%  within  the  first  year  for  a  diabetic  patient  with  a 
small  diameter  artery  undergoing  bare  metal  stent  deployment  (Leimgruber  et  ah,  1986). 
However,  there  is  currently  no  method  for  definitively  determining  if  restenosis  has 
occurred  without  an  invasive  procedure  to  determine  the  status  of  the  artery  in  question. 
While  these  procedures  have  become  low  risk  and  routine,  they  remain  fairly  expensive 
and,  as  with  all  procedures,  they  carry  some  small  but  real  risk  of  complication. 

One  possible  solution  for  this  issue  is  to  attach  a  MEMS  device  to  the  stent  placed 
into  the  artery.  This  device  would  contain  a  sensor  capable  of  taking  real  time 
measurements  of  either  blood  flow  or  pressure  and  transmitting  this  data  remotely  to  be 
analyzed.  By  monitoring  an  individual’s  blood  flow  through  the  artery,  doctors  would  be 
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able  to  notice  a  change  in  blood  flow  indicative  of  a  stenosis  prior  to  any  harmful  effects 
taking  place  and  without  any  of  the  risk  which  comes  from  an  additional  procedure. 

Creating  an  appropriate  device  would  require  four  parts:  a  sensor,  a  power  source, 
a  transmitter,  and  a  receiver.  Each  of  these  components  requires  specific  characteristics 
to  operate  within  the  artery  that  will  be  challenging  to  overcome.  The  first  of  these 
challenges  is  determining  a  sensor  capable  of  accurately  measuring  blood  flow.  The 
sensor  will  face  severe  size  limitations  due  to  the  relatively  small  diameter  of  coronary 
arteries  (<4.0  mm).  Not  only  will  it  need  to  be  small  enough  to  fit  on  the  stent,  but  it  also 
must  not  interfere  with  blood  flow  enough  to  function  as  a  blockage  of  its  own  or  pose 
any  risk  of  aggravating  pre-existing  coronary  conditions. 

Adding  a  traditional  battery-powered  source  will  require  more  space  to  be  utilized 
by  the  device  and  would  have  to  also  be  attached  to  the  stent,  further  increasing  the  risk 
of  the  device  interfering  with  blood  flow.  A  more  ideal  solution  would  be  to  have  the 
device  powered  by  the  flow  itself  (as  with  using  piezoelectric  material). 

As  with  the  other  components,  the  transmitter  will  need  to  be  attached  to  the  stent 
itself  in  such  a  manner  that  it  doesn’t  affect  blood  flow  or  increase  risk  for  clot  formation. 
If  the  device  were  to  transmit  constantly,  it  would  not  only  require  more  power,  but 
would  also  require  a  transmitter  strong  enough  to  be  capable  of  transmitting  through 
human  flesh,  which  creates  further  issues.  Another  possibility  would  be  a  passive  system 
that  is  queried  by  an  outside  source,  potentially  eliminating  the  need  for  an  active 
transmission  from  the  MEMS  device. 

A  passive  system  would  also  eliminate  certain  issues  with  the  receiver.  With  an 
active  transmitter  in  the  artery,  the  receiver  would  need  to  be  sensitive  enough  to  pick  up 
the  signal  of  the  transmitter,  which  would  be  weak  at  best.  With  a  query-type  system,  the 
majority  of  the  power  requirements  fall  to  the  receiver.  This  will  send  a  signal  to  the 
actual  device.  The  receiver  as  the  only  external  component  will  face  no  power 
restrictions,  so  placing  the  majority  of  the  power  requirements  on  this  component  can 
ease  the  design  requirements. 
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Due  to  the  complex  nature  of  MEMS  devices  and  the  stringent  requirements 
placed  on  the  device  by  the  coronary  circulation,  the  problem  is  an  ideal  candidate  for 
computer  simulations.  Simulations  allow  for  increased  flexibility,  enabling  the  analysis 
of  possible  components  in  various  locations  in  order  to  narrow  down  possibilities,  thus 
sparing  much  of  the  time  and  expense  that  goes  along  with  developing  and  carrying  out 
experiments. 

For  this  thesis,  COMSOL  multiphysics  software  (COMSOL  Inc.,  1998)  was 
selected  to  create  the  models  needed.  COMSOL  was  selected  for  its  ability  to  reproduce 
a  variety  of  physics  phenomena,  providing  flexibility  in  simulations  and  allowing  the 
system  and  device  to  be  modeled  and  evaluated.  Previous  COMSOL  simulations  of  flow 
through  stenoses  were  performed  by  Dean  (2013). 

Due  to  the  uncertainty  of  simulations,  models  will  be  compared  to  known 
experimental  data  in  order  to  validate  both  the  software  and  specific  models.  Once  the 
models  have  been  confirmed  as  accurate,  they  will  be  used  to  explore  different  options  to 
evaluate  how  both  the  system  and  MEMS  device  react  and  to  determine  the  best  course 
of  action.  Then  the  modified  system  will  be  created  for  experimentation  in  order  to  prove 
the  functionality  of  the  device  as  well  as  reconfirm  the  accuracy  of  the  model. 

B.  PRIOR  WORK 

1.  Experimental  Data 

The  initial  experimental  data  were  collected  by  Denardo  et  al.  (1994)  using  a 
scaled- up  model  of  coronary  blood  flow  (Figure  1)  and  a  Laser  Doppler  Velocimeter 
(LDV).  The  portion  of  the  model  that  was  being  tested  was  made  from  acrylic  tube  with 
an  internal  diameter  of  2.54  cm.  Clear  acrylic  was  selected  due  to  the  limitation  of  the 
LDV  that  requires  optical  clarity. 

In  order  to  ensure  laminar  flow  for  the  entrance  to  the  testing  region,  the  test 
section  length  was  233  cm.  A  constant  pressure  head  reservoir  with  an  overflow  pipe 
maintained  stable  pressure  and  continuous  flow.  A  gate  valve  was  adjusted  to  obtain  the 
desired  Reynolds  number.  Reynolds  number  was  calculated  as 
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Re=!^  (1) 

n 

o 

where  Re  is  the  unitless  Reynolds  number,  p  is  the  fluid  density  in  kg/m  ,  p  is  the 
dynamic  viscosity  in  Pa*s,  v  is  the  average  flow  velocity  and  D  is  the  inside  diameter  in 


meters. 


The  system  used  roughly  a  50%  mixture  of  glycerol  and  water  at  various  speeds 
in  order  to  maintain  constant  Reynolds  numbers  (Re)  at  200,  300  and  400.  Due  to  the 
changing  nature  of  the  mixture  in  use,  characteristics  of  the  fluid  were  determined 
through  the  use  of  standard  tables  prior  to  gathering  each  set  of  LDV  data. 

Acrylic  axisymmetric  stenoses,  modeled  in  Figure  2,  were  inserted  to  create 
constrictions  of  20,  40,  60,  75  and  90%  by  cross  sectional  area.  These  stenoses  began  at 
the  outside  wall  of  the  tube  and  angled  inward,  each  with  45  degree  edges  relative  to  the 
direction  of  the  tube  wall,  necessitated  in  order  to  gain  accurate  readings  from  the  LDV. 
Each  stenosis  was  two  cylinder  diameters  in  length  at  their  exterior,  regardless  of  their 
width. 


Constant 


Pressure  Head 

Reservoir  Overflow  Pipe 


Figure  1.  Scaled  model  of  artery  with  stenosis  used  for  experimental  data 

(after  Denardo,  1994). 
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Volume;  r-coordinate  (cm) 


Figure  2.  Schematic  representation  of  the  shape  of  the  90%  stenosis  with  a 

cutaway.  All  units  are  in  centimeters. 

Velocity  measurements  were  taken  at  nine  axial  locations.  For  this  thesis, 
negative  numbers  will  denote  upstream  locations,  and  the  axial  measurements  are  given 
in  terms  of  multiples  of  the  tube  diameter  relative  to  the  center  of  the  stenosis. 
Measurements  were  taken  upstream  at  locations  -1.500D,  -l.OOOD,  and  -0.475D,  at  the 
center  of  the  stenosis  (O.OOOD)  and  downstream  at  +0.475D,  +1.000D,  +2. GOOD,  and 
+2.500D.  At  each  of  these  positions,  measurements  were  taken  at  1  mm  radial 
increments.  The  measurements  at  the  +/-0.475D  marks  were  in  order  to  capture  the  flow 
at  the  edges  of  the  stenoses.  Focus  was  placed  on  downstream  flow  instead  of  upstream 
due  to  the  desire  to  capture  the  effects  of  the  blockage  downstream.  Cengel  and  Cimbala 
(2006)  provide  a  formula  for  the  axial  velocity  profile  u(r)  for  flow  in  a  pipe  given  by 

u(r)  =  2  X  Uavg  (l  -  y  (2) 
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where  is  the  average  fluid  velocity  in  m/s,  R  is  the  radius  of  the  cylinder  (1.27  cm 
for  the  experiment)  and  u(r)  is  the  velocity  as  a  function  of  the  radius  from  the  center  of 
the  tubing,  r.  The  quantity  u^vg  is  equal  to  half  the  maximum  velocity  for  fully  laminar 
flow.  This  equation  is  a  solution  to  the  Hagen-Poiseuille  equation  that  assumes  no-slip 
boundary  condition  at  the  pipe  wall  and  fully  developed  laminar  flow. 

All  of  the  data  points  were  taken  twice  for  each  stenosis,  once  without  and  then 
once  with  a  scaled  up  guidewire.  The  guidewire  in  use  was  a  0.5  cm  diameter  stainless 
steel  rod  inserted  into  the  center  of  the  tubing.  This  rod  was  painted  flat  black  so  as  to 
not  interfere  with  LDV  readings. 

The  LDV  used  to  collect  the  data  was  a  Lexus  Model  85  (Palo  Alto, 
CA)  with  a  single  color  argon  laser  using  a  wavelength  of  514.5  nm.  A  system  of  mirrors 
was  used  to  reflect  the  beam  into  the  laser  splitter  (TSI  Model  915  transmitter  optics,  TSI 
Model  918  receiving  optics).  The  photomultiplier  tube  (RCA  Model  4526,  Somerville, 
NJ;  DC  power  supply  Hewlett  Packard  Model  51  15 A,  Palo  Alto,  CA)  had  a  pinhole 
diameter  of  0.203  m.  A  TSI  Model  550  Intelligent  Flow  Analyzer  was  used  as  the  signal 
processor  to  analyze  the  data  in  real  time.  Net  radial  position  error  of  the  system  was 
within  ±  0.150  mm,  with  the  optics  placed  19.0  cm  from  the  sample  (Denardo  1994). 

2.  Modeling 

Building  off  the  work  of  Denardo  (1994),  Dean  (2013)  modeled  these  systems 
using  COMSOL  multiphysics  software.  Using  computer  models,  he  was  able  to  mimic 
each  of  the  10  possible  configurations  and  compare  the  results  to  the  experimental  200 
Reynolds  number  experimental  data  to  judge  the  accuracy  of  the  simulations. 

Due  to  the  axisymmetric  nature  of  the  system,  a  2D  model  was  used,  creating  only 
the  aspect  of  the  system  that  was  of  interest.  As  a  result,  the  model  consisted  of  only  a 
set  of  tubing  with  the  required  stenosis,  an  inlet,  an  outlet,  and  the  guidewire  (as 
necessary).  A  catheter  was  also  placed  within  the  model  for  studies  unrelated  to  his 
thesis.  (Dean,  2013)  For  this  reason,  the  catheter  (shown  in  Figure  3)  was  removed  in  all 
simulations  in  this  thesis.  In  order  to  ensure  parabolic  laminar  flow  occurred  before 
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reaching  the  stenosis,  the  required  entry  length  was  calculated  using  the  equation  given 
by  Cengel  and  Cimbala  (2006) 

Le^O.OSRexD  (3) 

where  Le  is  the  length  required  in  meters,  D  is  the  diameter  of  the  tube  in  meters,  and  Re 
is  the  Reynolds  number.  For  the  highest  Reynolds  number  of  400,  this  resulted  in  a 
requirement  of  61  centimeters  of  tubing  upstream  of  the  stenosis.  A  length  of  80 
centimeters  was  used  upstream  instead  in  order  to  allow  for  the  required  length  following 
the  interference  caused  by  the  catheter  (Dean  2013). 


Figure  3.  Model  inlet  with  catheter,  the  catheter  is  the  highlighted  section. 
The  stenosis  is  downstream  (not  pictured)  and  the  axis  of  symmetry  is 
located  at  r=0.  All  units  are  in  centimeters. 
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In  order  to  compare  the  model  to  the  experimental  data,  the  area  designated  as  the 
test  section  ranged  from  -1.5D  to  +2.5D  (Dean  2013)  encompassing  the  area  in  which 
measurements  were  recorded  by  Denardo  (1994).  A  zoom  in  of  this  area  is  shown  in 
Figure  4.  Due  to  radial  symmetry,  only  half  of  the  tube  needs  to  be  included  in  the 
model.  The  stenoses  were  modeled  after  those  in  the  in-vitro  experiment  and  share  the 
same  dimensions.  The  guidewire  in  these  models  was  assumed  to  be  both  perfectly 
straight  and  perfectly  centered  in  the  system  with  a  diameter  of  0.508  centimeters  (Dean, 
2013). 


Figure  4.  View  of  the  test  section  of  the  simulation  with  a  90%  stenosis  with 
guidewire.  The  axis  of  symmetry  is  located  at  r=0.  All  units  are  in 

centimeters. 

8 


Beyond  the  test  section,  an  additional  length  of  20  centimeters  of  tubing  was 
placed  in  the  model  to  show  the  downstream  effects  caused  by  the  stenosis.  At  the  end  of 
this  tubing,  an  outlet  boundary  condition  with  pressure  set  to  0  Pascals  gauge  pressure 
and  no  viscous  stress  was  imposed  on  the  system  in  order  to  simulate  an  indefinite  length. 

To  simulate  the  fluid  used  in  the  experiment,  a  custom  material  was  created  to 
allow  flexibility  in  matching  the  characteristics  of  the  mixture  used  in  the  experiment. 
The  Laminar  Flow  module  in  COMSOL  was  determined  to  be  adequate  to  model  the 
flow  that  occurred  in  the  experiment;  however,  certain  assumptions  were  required.  Most 
notably  the  fluid  was  assumed  to  be  incompressible,  and  the  boundary  conditions  initially 
were  all  assumed  to  be  no-slip  with  the  exception  of  at  the  guidewire  surface,  which  was 
given  a  moving  wall  boundary  condition.  A  moving  wall  boundary  condition  places  a  set 
velocity  on  the  fluid  at  the  point  of  contact  to  the  boundary.  It  achieves  this  velocity  at  the 
boundary  by  simulating  that  the  wall  has  infinite  length  and  it  moves  in  the  direction  of 
fluid  flow  with  a  constant  velocity  that  matches  the  desired  boundary  condition.  During 
the  course  of  modeling  the  system,  some  of  the  boundary  conditions  were  adjusted, 
which  will  be  addressed  later  in  this  thesis.  The  inlet  velocity  for  each  model  was 
derived  from  the  Reynolds  number  given  in  the  experiment. 


Dean’s  (2013)  simulations  at  Reynolds  number  200  showed  remarkably  close 
results  to  that  of  the  experimental  data  at  lower  grade  stenoses  (20  and  40%).  At  the 
higher  grade  stenoses  (60,  75  and  90%),  the  simulations  were  still  accurate,  although 
differences  were  more  apparent.  The  rapid  rise  of  velocity  at  the  entrance  to  the  stenosis 
created  some  issues  for  the  modeling  software,  particularly  in  the  cases  in  which  the 
guidewire  was  in  place.  The  previous  use  of  a  moving  wall  boundary  condition  for  the 
guidewire  also  provided  a  challenge  because  it  was  not  physically  realizable  and  proved 
to  be  less  effective  than  expected.  Thus,  a  slip  velocity  boundary  condition  was  used  at 
the  guidewire  surface  instead  for  all  stenoses.  The  formula  used  by  COMSOL  to  calculate 
the  value  of  slip  velocity  u  at  a  boundary  (assuming  zero  moving  wall  velocity  and  no 
thermal  creep  effects)  is  given  by: 


u  =  L. 


(4) 
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where  Lg  is  the  slip  length  in  meters  and  I— I  is  the  magnitude  of  the  derivative  of  the 

velocity  u  in  m/s  with  respect  to  radial  distance  r  in  meters  at  the  boundary.  COMSOL 
calculates  this  derivative  from  the  tangential  shear  stress  along  the  boundary  and  dividing 
it  by  the  fluid’s  viscosity  (COMSOL,  2013).  Figure  5  gives  a  visual  explanation  of  the 
slip  velocity  boundary  condition.  The  wall  on  the  left  (r  =  R)  has  a  no-slip  boundary 
condition  and  therefore  the  velocity  at  the  wall  goes  to  zero.  The  wall  on  the  right  (r  =  ro) 
has  a  partial  slip  boundary  condition.  By  assigning  a  specific  value  of  slip  length  Ls  and 
obtaining  the  value  of  the  slope  at  the  wall  from  the  COMSOL  numerical  solution  a  non¬ 
zero  velocity  at  the  wall  is  obtained. 


Figure  5.  Slip  velocity  boundary  condition. 

Dean  (2013)  showed  that  experimental  data  display  a  linear  relationship  between 
Reynolds  number  and  velocity  profiles  for  the  case  of  flow  through  a  stenosis  with  no 
guidewire  present.  Figure  6  re-creates  the  velocity  profile  obtained  by  Dean,  along  with 
the  velocity  profile  for  the  case  where  the  guidewire  is  present  for  the  75%  stenosis.  The 
linear  relationship  also  holds  true  for  the  guidewire  case. 
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Velocity  Profiles  of  75%  Stenosis  normalized  to  Re=100 


Velocity  Profiles  of  75%  Stenosis  normalized  to  Re=100 
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Figure  6.  Experimental  data  for  Re  200,  300  and  400  normalized  to  Re  100. 
The  left  profile  replicates  Dean’s  work.  The  right  profile  shows  that  the 
linear  relationship  in  velocity  also  holds  true  in  the  case  of  the 
guidewire  present. 


Dean  (2013)  was  unable  to  get  steady-state  results  for  the  90%  stenoses 
simulations  for  higher  Reynolds  numbers.  The  use  of  more  processing  power  and  the 
increase  in  the  number  of  iterations  required  for  the  COMSOL  solver  to  achieve  a  steady 
state  solution  enabled  the  simulations  to  run  for  both  Reynolds  number  300  and  400. 


C.  APPROACH 

This  thesis  focuses  on  the  next  step  in  developing  a  model  capable  of  emulating  a 
cardiovascular  system  by  exploring  the  necessary  features  of  COMSOL.  In  doing  so,  we 
will  evaluate  the  usefulness  of  COMSOL  as  a  modeling  software  to  understand  blood 
flow  in  constricted  arteries. 

This  thesis  will  also  address  the  weaknesses  in  the  previous  models  in  order  to 
continue  development  toward  the  end  goal  of  developing  a  MEMS  device  capable  of 
monitoring  blood  flow.  Most  notably,  focus  will  be  placed  on  the  higher- level  stenoses 
in  order  to  address  the  simulation  issues  previously  recognized.  An  emphasis  will  also  be 
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placed  on  addressing  the  boundary  condition  at  the  guidewire  surface  in  order  to  ensure 
accurate  modeling  in  all  systems. 
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II.  TIME  DEPENDENCY 


A.  IMPLEMENTATION 

Previously  completed  COMSOL  simulations  only  sought  steady  state  solutions. 
This  was  done  in  order  to  compare  the  simulations  to  the  experimental  data,  which  were 
taken  in  the  steady  state  condition.  In  order  to  continue  progress  toward  accurately 
modeling  a  human  cardiovascular  system  and  the  end  goal  of  developing  a  MEMs  device, 
a  model  capable  of  simulating  a  time-dependent  solution  is  necessary. 

Our  initial  attempt  at  time  dependency  involved  changing  none  of  the  previously 
set  variables  or  COMSOL  features  other  than  having  the  simulation  run  for  3  seconds 
while  capturing  a  result  every  0.1  seconds.  This  allows  the  ability  to  see  the  development 
of  the  system  prior  to  reaching  the  steady  state  solution.  Changing  only  one  variable  at  a 
time  also  allows  us  the  identification  of  any  new  effects  as  well  as  helping  to  troubleshoot 
should  any  issues  arise. 

Ultimately,  adding  time  dependency  created  new  issues  and  requirements  for 
COMSOL.  The  first  issue  encountered  was  the  inability  of  the  mesh  to  handle  time 
dependency  in  some  of  the  models.  COMSOL  uses  a  mesh  to  determine  specific  points 
at  which  to  calculate  a  solution,  and  interpolates  all  other  data  between  these  points. 
COMSOL  allows  for  specific  conditions  to  be  adjusted  manually  for  the  mesh,  however 
the  most  practical  and  efficient  method  for  our  models  was  to  use  the  preset  conditions. 

Using  the  preset  options  only  required  a  selection  of  how  fine  of  a  mesh  was 
needed  for  a  given  model.  Lor  these  models,  a  coarse  mesh  was  found  to  be  more  than 
capable  of  showing  the  necessary  definition.  Preset  options  allow  COMSOL  to  select 
where  to  place  the  points,  focusing  on  placing  more  points  in  areas  that  will  be  changing 
more  rapidly.  COMSOL  places  the  points  much  more  densely  near  the  stenosis, 
coinciding  with  the  area  of  interest  for  this  study,  as  seen  in  Ligure  7. 
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Figure  7.  Example  of  mesh  eonfiguration  in  a  90%  stenosis  with  guidewire 

model. 

Using  this  method  proved  very  sueeessful  for  steady  state  models;  however,  it 
ereated  eomplieations  when  time  dependeney  was  introdueed.  In  partieular,  for  the 
higher  grade  stenoses  (75%  and  90%)  with  guidewire,  the  model  was  unable  to  run  using 
the  previous  mesh  eonfiguration.  The  use  of  an  “adaptive  mesh”  feature  was  required  in 
order  to  run  these  time-dependent  solutions. 

The  adaptive  mesh  feature  of  COMSOL  is  only  available  for  time-dependent 

solutions.  The  user  selects  the  initial  mesh  setting  that  COMSOL  uses  for  the  first  time 

step.  From  that  point  on,  COMSOL  creates  a  new  mesh  for  each  given  time  step,  or  as 

specified  by  the  user.  While  this  adds  a  significant  amount  of  processing  time  to  the 

computation,  it  also  allows  COMSOL  to  shift  the  focus  of  the  computations  as  necessary 

to  capture  fluid  effects  as  they  develop.  This  is  of  particular  importance  for  any  fluid 
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effects  that  are  not  permanent  or  that  shift  throughout  the  system.  This  phenomenon  can 
be  seen  in  Figure  8,  where  the  adaptive  mesh  shifts  the  area  of  focus  for  computations  as 
an  effect  moves  downstream. 


Figure  8.  Example  of  self-adjustments  made  by  adaptive  mesh.  All  units  are 

in  centimeters. 
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B.  NEW  EFFECTS  IN  TIME-DEPENDENT  MODELS 

As  expected  for  the  lower  stenoses  (20,  40  and  60%),  the  effect  downstream  was 
fairly  insignificant.  For  the  higher  grade  stenoses,  the  effect  was  very  minimal  as  well 
with  no  guidewire.  This,  too,  was  expected,  as  the  reduction  in  area  was  far  less 
significant  than  the  75%  and  90%  with  guidewire  models. 

As  a  result,  the  primary  concern  with  this  effect  was  the  75%  and  90%  models 
when  the  guidewire  was  present.  Both  of  these  models  showed  fairly  significant  amounts 
of  disruption  in  the  initial  downstream  flow  following  the  stenosis,  although  this  was  far 
more  drastic  with  the  90%  (Figure  9). 

For  both  of  these  models,  however,  the  phenomenon  resulted  in  initial  eddies  that 
were  ultimately  pushed  through  the  system,  explaining  why  the  effect  was  never  seen 
during  the  steady  state  simulations. 
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Tlme-0,5  Surface:  Velocity  magr\ltude  (m/s) 


Time— 0.6  Surface:  Velocity  magnitude  (m/s) 


Tlme-0.7  Surface:  Velocity  magnitude  (m/s) 
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Figure  9.  Time  lapse  of  the  formation  and  movement  of  flow  using  a  steady 

input,  90%  stenosis,  and  guidewire. 

This  behavior  raised  more  questions.  Due  to  the  steady  state  nature  of  the  model 
in  previous  simulations,  the  disruption  downstream  was  fairly  insignificant.  It  was 
present  only  on  the  initialization  of  the  flow,  ultimately,  dissipating  in  time.  The  final 
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desired  system,  however,  requires  pulsatile  flow,  meaning  that  an  effect  on  an  increase  of 
flow  could  still  be  relevant  to  the  end  result,  and  the  ability  to  capture  this  effect  is 
essential  in  order  to  understand  how  the  system  behaves. 
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III.  COMPUTER  SIMULATIONS 


A.  20%  STENOSIS 

1.  No  Guidewire  Case 

In  order  to  achieve  a  20%  area  reduction  in  the  cylinder,  the  stenosis  was 
simulated  using  a  trapezoid  with  its  base  along  the  outer  wall  with  a  length  of  5.08  cm 
(two  tube  diameters).  Both  ends  leading  to  the  inner  wall  of  the  stenosis  had  an  angle  of 
45  degrees  to  achieve  a  height  of  0.13  cm  and  an  inner  length  of  4.69  cm  when  rounding 
of  the  corners  is  taken  into  consideration.  The  vertices  of  the  stenosis  wall  were  rounded 
by  means  of  a  fillet  geometry  function  with  a  radius  of  0.2  cm  to  approximate  the 
machining  characteristics  of  the  actual  acrylic  stenosis  used  in  the  experiment.  Rounding 
the  corners,  as  compared  to  using  sharp  corners,  had  no  effect  in  the  velocity  profiles  for 
any  of  the  stenoses  used.  However,  it  did  have  a  large  effect  in  the  shear  rate  magnitudes 
not  only  at  the  vertices  of  the  stenosis  but  throughout  the  stenosis  length  as  well.  A 
default  physics-controlled  mesh  was  used  with  a  “Fine”  element  size  (Dean  2013). 

For  the  Re300  (Reynolds  number  300)  20%  stenosis  with  no  guidewire  present, 
an  average  input  velocity  of  4.74  cm/s  was  chosen  (compared  to  an  average  input 
velocity  of  3.16  cm/s  obtained  in  the  Re200  experiment).  This  velocity  is  first 
approximated  by  setting  r  =  0  in  Eq.  (2)  to  find  the  average  velocity  of  the  fluid  from  the 
maximum  experimental  velocity  obtained  at  the  -1.5D  position.  The  average  velocity 
obtained  was  4.72  cm/s,  but  this  value  is  not  precise  due  to  the  errors  in  the  experimental 
measurement  of  the  velocity.  The  average  velocity  is  then  adjusted  to  better  match  the 
experimental  results.  The  dynamic  viscosity  for  the  experiment  was  4.4745  mPa*s,  and 
the  fluid  density  was  1 1 15.5  kg/m  .  A  no-slip  boundary  condition  is  enforced  at  all  points 
along  the  outer  wall  including  the  stenosis. 

Figure  10  shows  the  comparison  of  measured  and  simulated  velocity  profiles.  In 
order  to  plot  velocity  profiles  for  all  axial  locations  on  the  same  graph,  a  10  cm/s  velocity 
offset  was  added  to  each  progressive  profile  in  order  to  place  it  in  the  right  vertical  spot 
on  the  graph. 
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Figure  10.  Velocity  profile  for  the  20%  stenosis.  The  computer  simulation  is 
displayed  on  the  left  in  m/s.  Velocity  profiles  at  specific  axial  locations 
are  displayed  on  the  right  against  experimental  data  (blue  circles)  in 
cm/s.  Velocities  for  each  axial  location  are  spaced  by  10  cm/s  intervals. 


The  highest  velocity  measured  experimentally  was  10.15  cm/s  at  the  0  (mid¬ 
stenosis)  axial  location,  while  the  simulation  yielded  a  value  of  10.39  cm/s  for  the  same 
location.  A  good  agreement  between  the  measured  and  simulated  velocity  profiles  can  be 
seen  in  Figure  10. 

2.  Guidewire  Case 

The  setup  is  the  same  as  the  no  guidewire  case,  with  the  exception  of  the  inclusion 

of  the  guidewire  at  the  center  of  the  pipe.  Average  input  velocity  used  was  4.5  cm/s.  In 

this  case,  Eq.  (2)  will  provide  a  less  accurate  approximation  to  the  average  inlet  flow 

velocity  required,  given  that  the  measured  flow  at  the  guidewire  wall  does  not  go  to  zero. 
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The  same  phenomenon  of  non-zero  velocity  at  the  guidewire  wall  observed  by  Dean  in 
the  Re200  simulation  as  appears  in  the  ReSOO  simulation,  albeit  at  a  higher  magnitude. 


Figure  11.  Velocity  profile  for  the  20%  stenosis  with  guidewire  present. 

Computer  simulation  is  displayed  on  the  left  is  m/s.  Velocity  profiles  in 
cm/s  displayed  on  the  right  against  experimental  data  (blue  circles). 


The  experimental  results  in  Figure  11  show  that  the  velocity  of  the  flow  at  the 
point  closest  to  the  guidewire  at  the  -1.5D  location  is  1.55  cm/s,  at  -1.0  the  velocity  is 
1.41cm/s,  and  inside  the  stenosis  the  velocity  closest  to  the  guidewire  surface  goes  up  to 
3.63,  3.01  and  2.6  cm/s  for  the  -.475D,  O.OD,  and  -I-.475D  locations,  respectively.  In  order 
to  simulate  this  non-zero  velocity  phenomena  observed  at  the  guidewire  surface,  a  slip- 
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velocity  boundary  condition  was  set  along  the  guidewire  wall  using  a  slip  length  Ls  of  1.7 
mm.  The  maximum  experimental  velocity  recorded  for  this  stenosis  was  8.86  cm/s  at  the 
O.OD  axial  location,  compared  with  a  value  of  8.26  cm/s  obtained  from  the  simulation  at 
the  same  axial  location.  A  good  agreement  between  the  measured  and  simulated  velocity 
profiles  can  be  seen  in  Figure  1 1 . 

B.  40%  STENOSIS 

1.  No  Guidewire  Case 

The  dynamic  viscosity  and  density  of  the  fluid  used  for  the  40%  stenosis 
simulation  were  the  same  values  as  the  20%  stenosis  given  that  the  same  fluid  was  used 
for  the  experiment.  The  40%  constriction  is  achieved  by  giving  the  stenosis  a  height  of 
0.29  cm.  Since  the  outer  length  of  the  stenosis  and  the  angles  are  fixed  to  two  diameter 
lengths  and  45  degrees  respectively,  the  inner  length  of  the  stenosis  is  4.21  cm.  Fillets  of 
radius  2  cm  were  used  to  round  the  corners  of  the  vertices.  Figure  12  displays  the 
measured  and  simulated  velocity  profiles. 
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Figure  12.  Velocity  profile  for  the  40%  stenosis  with  no  guidewire.  Velocities 
for  each  axial  location  are  spaced  by  10  cm/s  intervals. 


The  maximum  experimental  velocity  recorded  was  12.32  cm/s  at  the  +.475D  axial 
location  compared  to  a  value  of  12.23  cm/s  obtained  from  the  simulation  at  the  same 
location.  A  good  agreement  between  the  measured  and  simulated  velocity  profiles  can  be 
seen  in  Figure  12. 

2.  Guidewire  Case 

The  fluid  properties  were  the  same  as  the  ones  used  for  the  no-guidewire  case.  A 
constant  moving  wall  velocity  at  the  guidewire  was  not  able  to  replicate  the  experimental 
results  given  that  a  larger  velocity  was  required  near  the  guidewire  wall  at  the  stenosis 
height.  However,  the  slip  velocity  boundary  condition  with  viscous  slip  with  slip  length 
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of  0.88  mm  was  able  to  replicate  the  experimental  results  at  the  guidewire.  Figure  13 
shows  the  comparison  of  measured  and  simulated  velocity  profiles.  There  is  still  a 
discrepancy  at  the  -.475D  position  close  to  the  stenosis  wall  that  is  due  to  the  no-slip 
boundary  condition  forced  at  the  stenosis  wall.  The  40%  stenosis  with  guidewire 
simulation  is  the  smallest  stenosis  to  form  a  recirculation  region  past  the  exit  of  the 
stenosis,  where  a  maximum  negative  velocity  of  -3.27  cm/s  occurs. 
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Figure  13.  Velocity  profile  for  the  40%  stenosis  with  guidewire.  Experimental 
data  show  a  slip  boundary  condition  at  the  entrance  of  the  stenosis, 
while  the  no-slip  boundary  condition  imposed  on  the  simulation  forces 

the  velocity  to  go  to  zero. 

The  40%  stenosis  forces  the  experimental  velocity  recorded  closest  to  the 

guidewire  wall  to  increase  from  4.03  cm/s  at  the  -l.OD  location  to  6.25,  6.92,  and  7.29 
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cm/s  for  the  -.475D,  O.OD  and  +.475D  locations,  respectively.  The  slip  velocity  boundary 
condition  established  at  the  guidewire  surface  is  able  to  mimic  this  rise  in  velocity  inside 
the  stenosis,  where  clearly  a  constant  moving  wall  velocity  would  fall  short  of 
accomplishing  the  same  result.  The  maximum  experimental  velocity  observed  was  11.72 
cm/s  at  the  +.475D  location,  midway  between  the  stenosis  and  guidewire  walls.  The 
simulation  velocity  obtained  at  the  same  point  was  1 1.07  cm/s. 

C.  60%  STENOSIS 

1.  No  Guidewire  Case 

The  60%  constriction  is  achieved  with  a  stenosis  height  of  0.47  cm,  which  yielded 
an  inner  stenosis  wall  length  of  4.14  cm.  The  dynamic  viscosity  for  the  experiment  was 
6.2536  niPa*s  and  fluid  density  of  1133.2  kg/m  .  A  no-slip  boundary  condition  is 
enforced  at  all  points  along  the  outer  wall  including  the  stenosis.  Input  velocity  had  a 
value  of  6.1  cm/s.  The  60%  stenosis  is  the  lowest  grade  stenosis  to  yield  a  negative 
velocity  due  to  backflow  when  no  guidewire  is  present.  Figure  14  shows  the  measured 
and  simulated  velocity  profiles. 
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Figure  14.  Velocity  profile  for  the  60%  stenosis.  Experimental  values  show  a 
non-zero  velocity  at  the  surface  of  the  stenosis.  Velocities  for  each  axial 
location  are  spaced  by  20  cm/s  intervals. 


Backflow  occurred  past  the  stenosis  and  reached  a  maximum  speed  of  1.12  cm/s. 
The  highest  fluid  velocity  from  the  simulation  was  23.73  cm/s  and  occurred  centerline 
between  the  -I-.475D  and  -i-l.OD  locations,  coinciding  with  the  location  where  the  stenosis 
ended.  The  maximum  value  obtained  from  the  experiment  occurred  at  the  -i-l  .OD  location 
closest  to  centerline  with  a  value  of  23.08  cm/s,  versus  23.64  cm/s  obtained  by  the 
simulation  at  the  same  location. 
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2. 


Guidewire  Case 


The  same  fluid  from  the  no  guidewire  case  was  used.  The  average  input  velocity 
was  6.5  cm/s  and  a  viscous  slip  boundary  condition  with  slip  length  1.5  mm  was  used  at 
the  guidewire.  An  exponentially  decaying  moving  wall  boundary  condition  was 
attempted  in  order  to  replicate  the  behavior  near  the  stenosis  wall  at  the  -.475D  position. 
This  exponentially  decaying  velocity  function  was  successful  at  replicating  the 
experimental  results.  However,  the  experimental  data  does  not  include  values  for  the 
velocity  close  to  the  stenosis  wall.  As  a  result,  many  different  values  for  the  amplitude 
and  decay  constant  for  the  exponentially  decaying  velocity  function  could  replicate  the 
experimental  results.  A  no-slip  boundary  condition  was  decided  for  all  stenoses  until  a 
new  experiment  is  performed  with  better  velocity  resolution  near  the  stenosis  walls. 
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Figure  15.  Velocity  profile  for  the  60%  stenosis  with  guidewire. 


The  addition  of  the  guidewire  increased  the  maximum  fluid  velocity  from  23.73 
cm/s  to  26.83  cm/s  at  a  location  0.4  cm  away  from  centerline  and  a  position  of  +0.7D. 
The  highest  velocity  measured  experimentally  occurred  at  the  +.475D  position  with  a 
value  of  26.52  cm/s,  compared  with  26.06  cm/s  from  the  simulation  at  the  same  location. 
The  highest  velocity  from  the  simulation  occurred  at  the  +.7D  location  with  a  value  of 
26.83  cm/s. 
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D.  75%  STENOSIS 


1.  No  Guidewire  Case 

The  75%  stenosis  constriction  is  achieved  with  a  stenosis  height  of  0.635  cm, 
which  yielded  an  inner  stenosis  wall  length  of  3.81  cm.  The  dynamic  viscosity  for  the 
experiment  was  11.833  mPa*s  and  fluid  density  of  1167  kg/m  .  Input  velocity  for  the 
simulation  had  a  value  of  13.8  cm/s.  Due  to  the  small  cross-sectional  area  caused  by  the 
presence  of  a  large  stenosis  the  number  of  measurements  inside  the  stenosis  was  doubled, 
as  shown  in  Figure  16.  This  was  done  by  measuring  at  a  distance  of  half  a  millimeter 
between  measurements  compared  to  one  millimeter  for  lower  stenoses. 


Figure  16.  Velocity  profile  for  the  75%  stenosis  with  no  guidewire. 
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The  reason  the  velocity  profile  throughout  the  channel  is  higher  in  the  simulation 
for  the  -.475D  position  is  because  of  the  no-slip  boundary  condition  imposed  at  the 
stenosis  wall,  which  by  conservation  of  mass  forces  the  flow  to  have  a  higher  velocity 
throughout  the  rest  of  the  channel.  The  highest  simulation  velocity  had  a  value  of  79  cm/s 
at  the  -I-0.7D  diameter  position.  The  maximum  experimental  velocity  of  78.81  cm/s 
occurred  at  the  -I-.475D  position.  The  velocity  obtained  from  the  simulation  for  the  same 
point  was  78.41  cm/s. 

2.  Guidewire  Case 

In  contrast  to  all  previous  stenosis  data,  the  75%  stenosis  experimental  data 

required  a  normalization  factor  to  account  for  the  fact  that  three  different  fluids  with 

different  densities  and  viscosities  were  used  for  the  experiment.  Fluid  one  had  a  viscosity 

-2 

of  11.833  mPa*s,  a  density  of  1167  kg/m  ,and  was  used  for  the  -1.5D  measurement. 

-2 

Fluid  two  had  a  viscosity  of  4.175  mPa*s,  density  of  1112  kg/m  and  was  used  for  the  - 
-1-1. OD  measurement.  Fluid  three  had  a  viscosity  of  5.119  mPa*s,  a  fluid  density  of  1123 
kg/m^  and  was  used  for  the  -.475D,  O.OD,  -I-.475D,  -i-l.OD  and  -I-2.5D  measurements. 
Previous  attempts  at  matching  the  simulation  to  experimental  data  failed  because  only 
fluid  three  was  used  for  the  simulation.  It  was  discovered  that  although  different  fluids 
may  be  normalized  by  means  of  taking  the  ratio  of  their  viscosities  and  densities,  the  slip 
velocity  at  the  guidewire  wall  depends  non-linearly  on  the  fluid’s  viscosity.  Therefore  a 
different  simulation  was  required  for  each  fluid.  The  simulation  for  fluid  one  used  an 
input  velocity  of  5.67  cm/s  and  a  viscous  slip  boundary  condition  at  the  guidewire  surface 
with  a  slip  length  of  2.5  mm;  the  simulation  for  fluid  two  used  an  input  velocity  of  4.58 
cm/s  and  a  slip  length  of  2.5  mm,  and  the  simulation  for  fluid  three  used  an  input  velocity 
of  5.72  cm/s  and  a  slip  length  of  .5  mm.  To  be  able  to  plot  the  experimental  data  for  all 
locations  on  the  same  graph,  the  results  for  the  -1.5D  and  -l.OD  values  were  normalized 
to  fluid  three’s  values.  Figure  17  shows  the  comparison  of  measured  and  simulated 
velocity  profiles. 
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Figure  17.  Velocity  profile  for  the  75%  stenosis  with  guidewire. 


The  highest  experimental  velocity  obtained  was  39.43  cm/s  at  the  +.475D  location 
compared  with  a  simulation  value  of  35.15  cm/s  at  the  same  location.  The  largest 
discrepancy  occurs  at  the  stenosis  surface  for  the  -.475D  location,  and  it  is  due  to  the  no¬ 
slip  boundary  condition  imposed  on  the  stenosis  surface  which  forces  the  velocity  to  go 
to  zero. 


31 


E.  90%  STENOSIS 


1.  No  Guidewire  Case 

The  90%  stenosis  constriction  is  achieved  with  a  stenosis  height  of  0.86  cm, 
which  yielded  an  inner  stenosis  wall  length  of  3.36  cm.  The  90%  stenosis  without 
guidewire  experimental  setup  required  two  fluids.  Fluid  one  was  used  for  the  -1.5D,  - 
l.OD  and  -.475D  positions,  with  a  viscosity  value  of  11.833  mPa*s  and  density  of  1167 
kg/m^.  Fluid  two  was  used  for  the  O.OD,  +.475D,  +1.0D,  +1.5D,  +2.0D  and  +2.5D 
positions  with  a  viscosity  of  4.175  mPa*s  and  density  of  1112  kg/m  .  Normalization  of 
experimental  values  was  required  as  was  the  case  for  the  75%  stenosis  with  guidewire. 
The  simulation  was  run  using  fluid  two’s  properties.  The  average  input  velocity  was  4.5 
cm/s.  Figure  18  shows  the  measured  and  simulated  velocity  profiles. 


0  1^0 


Figure  18.  Velocity  profile  for  the  90%  stenosis. 
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For  the  90%  stenosis  without  guidewire  setup,  the  maximum  velocity  did  not 
occur  at  a  location  that  matched  the  location  of  the  stenosis.  Instead,  the  maximum 
experimental  velocity  occurred  at  the  +1.5D  axial  location,  with  a  value  of  60.49  cm/s, 
compared  with  61.78  cm/s  obtained  from  the  simulation  at  the  same  location. 

2.  Guidewire  Case 

For  the  90%  stenosis  with  guidewire,  two  different  fluids  were  used  and 
normalization  was  performed  to  graph  the  results.  Fluid  one  had  a  viscosity  of  4.4745 
mPa*s  and  density  of  1 1 10  kg/m  .  It  was  used  for  the  experimental  measurements  at  the 
-1.5D,  -l.OD  and  -.475D  positions.  Fluid  two  had  a  viscosity  of  4.175  mPa*s  and  a 
density  of  1112  kg/m  .  Similar  to  the  75%  stenosis  with  guidewire  case,  a  simulation  had 
to  be  run  for  each  fluid  to  account  for  the  non-linear  characteristics  of  the  slip  velocity  at 
the  guidewire  surface.  Fluid  one  used  an  input  velocity  of  4.7  cm/s  and  a  slip  velocity 
boundary  condition  at  the  guidewire  wall  with  slip  length  of  0. 1  mm  and  fluid  two  used 
an  input  velocity  of  5.2  cm/s  with  a  slip  length  of  0.3  mm.  The  outer  wall  (including  the 
stenosis)  had  a  no-slip  boundary  condition.  From  Figure  19  it  is  clear  that  the  simulation 
cannot  match  the  experimental  results  at  the  guidewire  wall. 
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Figure  19.  Velocity  profile  for  the  90%  stenosis  with  guidewire. 


The  velocity  of  the  fluid  near  the  guidewire  in  the  90%  simulation  with  guidewire 
is  smaller  than  the  experimental  results.  Contrary  to  the  75%  simulation  with  guidewire, 
the  90%  simulation  is  not  able  to  match  experimental  results  properly  for  the  downstream 
axial  locations.  This  is  caused  by  the  no-slip  boundary  condition  which  forces  the 
velocity  at  the  stenosis  wall  to  be  zero,  therefore  increasing  the  flow  near  the  guidewire. 
By  the  time  the  fluid  exits  the  stenosis,  the  fluid  flux  is  still  concentrated  in  the  central 
regions  of  the  channel,  so  the  magnitude  of  the  velocity  has  to  be  higher  in  the  open  area 
of  the  channel,  overshooting  the  experimental  results.  The  75%  simulation  is  still  able  to 
match  experimental  results  downstream  because  by  the  time  the  fluid  passes  the  -I-.475D 
location  the  experimental  fluid  velocity  near  the  stenosis  wall  approaches  zero  and 

matches  the  no-slip  boundary  condition  of  the  simulation.  The  same  shortcoming  was 

34 


observed  by  Dean  in  the  Re200  simulation  and  again  by  the  authors  when  performing  the 
simulation  for  Re400  values.  The  maximum  experimental  velocity  observed  was  104.23 
cm/s  at  the  mid-stenosis  axial  location,  compared  with  a  simulation  velocity  of  101.77 
cm/s  at  the  same  location. 

F.  SHEAR  RATE 

An  analysis  of  shear  rate  was  performed  using  data  from  the  stenosis  simulations 
for  both  the  guidewire  surface  (when  applicable)  and  the  outer  wall  surface  for  all 
stenoses.  Shear  rate  is  important  in  cardiovascular  medicine  because  it  affects 
atherosclerosis  (Giddens,  Zarins  &Glagov,  1993),  thrombosis  (Badimon,  Chesebro  & 
Badimon,  1992),  and  restenosis  after  coronary  intervention  (Koskinas,  Chatzizisis, 
Antoniadis  &  Giannoglou,  2012).  A  distance  of  1  micron  from  the  surface  was  chosen  for 
all  the  shear  rate  measurements.  Shear  rate  values  were  readily  obtained  from  the 
simulations.  COMSOL  analyzes  the  shear  stress  from  the  numerical  simulations  and  the 
fluid’s  viscosity  to  obtain  the  shear  rate  at  the  desired  position.  Unfortunately,  the  results 
obtained  cannot  be  compared  against  the  experiment  data  due  to  the  physical  limitations 
of  the  measuring  device.  The  LDV  was  not  able  to  obtain  precise  measurements  close  to 
the  surfaces,  so  meaningful  shear  rate  data  cannot  be  extracted  from  the  experiment.  The 
guidewire  surface  had  a  slip  velocity  boundary  condition,  with  a  slip  length  magnitude 
that  matched  the  conditions  already  stated  for  each  simulation. 

1.  Guidewire  Surface 

It  is  clear  from  Figure  20  that  the  magnitude  of  the  shear  rate  at  the  guidewire 
surface  increases  as  the  stenosis  grade  increases,  not  only  at  the  axial  location  where  the 
stenosis  is  located  (-ID  to  -I-ID)  but  throughout  the  full  length  of  the  guidewire.  The 
maximum  shear  rate  magnitude  obtained  matched  the  section  of  the  guidewire  at  an  axial 
distance  coinciding  with  the  upstream  vertex  of  the  stenosis.  The  shear  rate  at  the 
guidewire  for  20%  stenosis  reached  a  maximum  value  of  29.73  s'\  which  is  a  27% 
increase  from  the  nominal  shear  rate  value  downstream  of  the  stenosis.  The  40%  stenosis 
had  a  maximum  shear  rate  magnitude  of  57.28  s’\  which  was  twice  as  high  as  the 
nominal  downstream  shear  rate  value  of  27.75  s’\  The  60%  stenosis  had  a  maximum 
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shear  rate  magnitude  value  of  113.7s'\  Similarly  to  the  75  and  90%  stenoses  there  is  a 
gradual  decrease  in  shear  rate  downstream  of  the  stenosis,  but  the  simulation  setup  is  not 
long  enough  to  achieve  a  constant  shear  rate  magnitude.  The  75%  stenosis  simulation 
yielded  a  maximum  shear  rate  at  the  guidewire  of  258.85  s'^  and  the  90%  maximum  was 
3200  s'\  The  color  scale  in  Figure  20  was  chosen  so  that  any  value  of  shear  rate  above 
200  s'^  would  show  as  a  single  color  in  order  to  maximize  the  contrast  of  the  lower 
magnitude  values  so  the  phenomena  downstream  of  the  guidewire  would  be  visible. 


R=300  Shear  rate  along  guidewire  (slip  velocity  b/c) 


Axial  Location  (Diameters) 


Figure  20.  Top:  Logarithmic  graph  of  shear  rate  magnitude  1  micron  away  from 
the  guidewire  units  in  1/s.  Shown  in  red  is  the  location  of  the  stenosis 
along  the  axial  location.  Bottom:  Shear  rate  magnitude  profile  for  the 
90%  stenosis  with  guidewire. 
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Figure  20  illustrates  that  the  90%  stenosis  simulation  exhibits  fluctuation  of  shear 
rate  past  the  +4D  mark  downstream  of  the  stenosis.  The  velocity  profile  of  the  fluid 
shows  a  recirculation  zone  that  bends  the  flow  towards  the  outer  wall  starting  at  the 
+4.5D  location.  As  the  flow  shifts  away  from  the  guidewire,  the  shear  rate  at  the 
guidewire  surface  decreases  dramatically  since  the  slope  of  the  velocity  is  less  steep  due 
to  a  recirculation  region  that  emerges  at  the  guidewire  at  the  +5D  location,  and  increases 
in  the  space  between  the  guidewire  and  outer  wall  as  the  velocity  increases  in  magnitude 
and  its  slope  is  steeper. 

The  same  shear  rate  analysis  was  performed  using  Dean’s  moving  wall  boundary 
condition  at  the  guidewire  for  comparison  purposes.  Even  for  the  smallest  of  the  stenoses 
the  maximum  shear  rate  magnitude  varied  by  as  much  as  22%.  Though  either  boundary 
condition  was  adequate  to  replicate  experimental  results  of  velocity  profiles,  the  moving 
wall  boundary  condition  provides  inadequate  results  in  shear  rate.  The  maximum  shear 
rate  recorded  at  the  guidewire  was  3200  s'^  which  occurred  for  the  90%  stenosis  case 
using  a  slip  velocity  boundary  condition  at  the  guidewire  surface,  compared  to  4661  s’^ 
for  the  moving  wall  boundary  condition. 

2.  Outer  Wall  with  No  Guidewire 

Figure  21  shows  the  results  of  measuring  the  magnitude  of  the  shear  rate  one 
micron  away  from  the  outer  wall  for  the  simulations  with  no  guidewire  present.  There  is  a 
clear  pattern  inside  the  stenosis  wall  of  increased  shear  rate  as  the  size  of  the  stenosis  is 
increased.  The  maximum  value  of  the  shear  rate  always  occurred  at  the  upstream  vertex 
of  the  stenoses,  as  shown  in  Figure  21. 
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R=300  Shear  rate  along  outer  wall  (no  guidewire,  no  slip  b/c  at  outer  wall) 


Figure  21.  Logarithmic  graph  of  shear  rate  along  the  outer  wall  for  simulations 
with  no  guidewire  present.  No  slip  boundary  condition  along  the  length 
of  the  outer  wall  (including  stenosis).  Units  in  1/s. 


The  maximum  shear  rates  obtained  were  63.68,  141.10,  749.45,  1433.04  and 
1627.85  s'^  for  the  20,  40,  60,  75  and  90%  stenoses,  respectively.  Downstream  from  the 
stenoses  the  pattern  is  more  complicated.  The  40%  stenosis  simulation  develops  a  dip  in 
shear  rate  at  the  +1.7D  location,  while  the  60%  stenosis  simulation  develops  a  similar  dip 
at  the  +4.4D  location.  The  simulation  ended  at  100  cm  to  mimic  the  experimental  setup, 
so  the  dips  in  shear  rate  for  the  75%  and  90%  stenoses  are  not  visible  in  the  figure; 
however,  subsequent  simulations  proved  that  the  dips  do  occur  further  downstream.  For 
the  90%  stenosis  data,  the  shear  rate  just  downstream  of  the  stenosis  has  a  smaller 
magnitude  than  the  rest  of  the  stenoses,  indicating  a  region  of  near  zero  flow  where  the 
gradient  of  the  velocity  is  fairly  constant  near  the  wall.  The  maximum  shear  rate 
downstream  of  the  stenosis  occurred  for  the  75%  stenosis  simulation,  with  a  maximum 
shear  rate  of  52.25  s'^  at  the  +3. 3D  location.  The  maximum  downstream  value  for  60% 
stenosis  was  14.86  s'^  at  +1.6D,  and  for  the  90%  stenosis  the  maximum  shear  rate  value 
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obtained  was  29.62  s'^  at  the  +6. 3D  location.  Maxima  for  the  20  and  40%  stenoses 
occurred  at  points  further  out  from  the  simulation  limit. 


3.  Outer  Wall  with  Guidewire 

Figure  22  shows  the  results  of  measuring  the  magnitude  of  the  shear  rate  one 
micron  away  from  the  outer  wall  for  the  simulations  with  the  guidewire  present.  The 
presence  of  the  guidewire  increases  the  shear  rate  at  the  stenosis  wall  for  all  simulations, 
as  shown  in  Figure  22. 


R=300  Shear  rate  along  outer  wall  (with  guidewire, 
no  slip  B/C  at  outer  wall,  slip  velocity  b/c  at  guidewire) 


Axial  Location  (Diameters) 


Figure  22.  Logarithmic  graph  of  shear  rate  along  the  outer  wall  for  simulations 
with  guidewire  present.  No  slip  boundary  condition  along  the  length  of 
the  outer  wall  (including  stenosis).  Units  in  1/s. 


The  maximum  value  of  the  shear  rate  obtained  was  82.63,  169.80,  479.12,  706.32 
and  4315.78  s'^  for  the  20,  40,  60,  75  and  90%  stenoses,  respectively.  The  maximum 
value  always  occurred  at  the  upstream  vertex  of  the  stenoses.  Compared  to  the  previous 
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study,  the  presence  of  the  guidewire  causes  the  75%  downstream  dip  in  shear  rate  to 
travel  upstream  to  the  +6. 2D  location,  compared  to  +10.4D  for  the  no  guidewire  case. 
The  90%  stenosis  simulation  shows  a  discrepancy  that  matched  the  behavior  seen  at  the 
same  axial  location  in  the  guidewire  as  seen  in  Figure  20.  Fluid  flow  converging  from  the 
middle  of  the  channel  strikes  the  outer  wall  at  the  +5. 2D  location,  causing  areas  of  high 
shear  rate  at  the  outer  wall  that  reach  magnitudes  of  up  to  975.22  s’\ 
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IV.  PULSATILE  FLOW 


A.  SIMULATING  A  HEARTBEAT 

This  thesis  offers  only  an  initial  investigation  of  pulsatile  flow  mimicking  the 
human  heartbeat.  For  further  development,  the  velocity  profile  required  to  model  a 
human  heartbeat  must  be  established.  After  an  extensive  literature  search,  it  was 
determined  that  a  significant  amount  of  future  work  will  be  required  to  determine  the 
various  velocity  profiles  necessary  for  modeling. 

Not  only  are  there  variations  in  the  profiles  from  person  to  person,  but  the 
velocity  profile  can  also  change  drastically  depending  on  the  location  within  the 
cardiovascular  system  being  monitored.  The  first  of  these  two  issues  will  require  further 
research  in  determining  a  standard  example  to  use  for  this  system.  Ultimately,  however, 
if  a  sensor  is  going  to  be  developed  for  medical  use,  it  must  be  functional  for  the  full 
range  of  profiles  that  may  be  encountered. 

In  order  to  address  the  second  issue  of  location  within  the  cardiovascular  system, 
it  will  be  essential  to  determine  velocity  profiles  in  the  areas  applicable  to  these  medical 
procedures.  While  the  nature  of  this  problem  results  in  a  limited  portion  of  the 
cardiovascular  system,  further  research  will  have  to  be  done  specifying  these  locations, 
determining  the  individual  profiles,  and  ensuring  the  sensor  is  capable  of  functioning  in 
various  environments. 

For  the  purpose  of  this  thesis,  we  sought  only  a  proof  of  concept  for  pulsatile  flow 
to  determine  what  issues  are  inherent  with  pulsatile  flow  being  introduced  to  previous 
models.  To  complete  this  aspect  of  the  model,  we  required  a  variable  inlet  velocity  in 
the  general  shape  of  a  heartbeat  to  simulate  pulsatile  flow. 

Using  COMSOL,  we  are  able  to  enter  as  many  or  few  points  as  necessary  with  no 
requirement  for  the  points  to  be  equidistant  in  the  time  scale.  COMSOL  then  uses  a 
piecewise  cubic  solution  to  interpolate  between  the  points.  The  result  is  the  function  seen 
in  Figure  23,  named  beat  in  the  models.  This  function  was  used  only  as  a  method  to  scale 
previous  inputs  and  does  not  represent  the  actual  inlet  velocity  in  meters  per  second. 
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Beat(t) 


Beat(t) 


Figure  23.  Unitless  funetion  applied  to  pulsatile  flow  to  mimic  a  heartbeat. 

Time  is  in  seconds. 

Pulsatile  flow  also  requires  the  use  of  COMSOL’s  mod  function.  This  allows  for 
a  function  (such  as  the  beat  function)  to  be  repeated,  for  the  purposes  of  a  these 
simulations  it  was  repeated  once  every  second.  This  is  critical  in  order  to  establish  how 
the  system  will  function  in  a  realistic  situation  over  extended  periods  of  time.  By 
repeating  the  pulse  we  are  able  to  determine  the  effects  of  a  continuous  system. 
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B.  RESULTS 


As  was  expected,  following  the  original  time-dependent  simulations  the 
constriction  caused  by  the  higher  stenoses  when  the  guidewire  was  present  did  result  in 
eddies  forming,  as  seen  in  Figure  24.  These  eddies  virtually  dissipated  between  pulses 
and  did  not  continue  to  flow  downstream,  but  they  are  significant  in  the  area  immediately 
following  the  stenosis. 


Figure  24.  Time  lapse  of  downstream  effect  in  a  90%  stenosis  model  with 

guidewire  using  pulsatile  flow.  Each  image  represents  a  0.1  s  time  step 

during  a  pulse. 
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It  was  found  that  the  presence  of  the  guidewire  played  an  even  larger  role  in  the 
formation  of  these  eddies  than  the  size  of  the  stenosis.  This  is  because  of  the  effect  it 
plays  in  changing  the  opening  area  of  the  tubing.  While  the  reduction  in  stenosis  has  a 
greater  effect  on  the  overall  area,  the  guidewire  separates  the  remaining  area  into  two 
parts,  forcing  a  more  drastic  affect. 

The  result  is  that  without  a  guidewire,  even  the  90%  stenosis  shows  virtually  no 
effect  downstream  as  seen  in  Figure  25.  Further  simulations  will  have  to  be  completed 
with  a  more  accurate  model  of  a  heartbeat,  but  these  results  show  that  it  is  likely  that  a 
higher  level  stenosis  may  have  this  effect  only  when  the  guidewire  is  inserted.  Without 
the  guidewire  implanted  (as  will  be  the  normal  case),  the  effect  is  most  likely  going  to  be 
negligible. 


Figure  25.  Time  lapse  of  downstream  flow  in  a  90%  stenosis  model  without 
guidewire  using  pulsatile  flow.  Each  step  requires  0.1  second  time  step 

in  a  pulse. 
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V.  CONCLUSIONS  AND  FUTURE  WORK 


Higher  processing  power  and  an  increase  in  the  number  of  iterations  required  by 
the  software  solver  enabled  COMSOL  to  achieve  a  steady  state  solution  for  all  three 
Reynolds  numbers  (200,  300  and  400).  A  non-zero  velocity  at  the  guidewire  wall  caused 
by  a  partial  slip  condition  was  verified  in  the  experimental  results  for  all  Reynolds 
numbers.  The  use  of  a  “slip  velocity”  boundary  condition  was  proven  to  be  successful  at 
simulating  experimental  results  at  the  guidewire  surface  for  all  stenosis  grade  conditions 
without  the  need  for  the  artificial  “moving  wall”  boundary  condition  used  in  previous 
studies,  which  required  arbitrary  velocity  magnitudes  at  different  points  along  the 
guidewire  to  match  experimental  results. 

The  shear  rate  study  at  the  boundaries  provided  insights  into  the  shear  stresses 
encountered  by  the  guidewire  and  outer  walls  and  how  their  magnitude  changes  as  a 
function  of  level  of  constriction.  By  running  simulations  with  both  moving-wall  and  slip- 
velocity  boundary  conditions  at  the  guidewire  wall,  it  was  possible  to  see  the  difference 
in  shear  rate  magnitude  encountered  when  using  each  boundary  condition.  Possible 
causes  for  the  non-zero  fluid  velocity  at  the  guidewire  include  the  formation  of  a 
boundary  layer  due  to  the  non-hydrophobic  paint  used  for  the  guidewire  or  a  partial  slip 
condition  that  occurs  due  to  the  irregularities  in  the  guidewire  surface.  A  new  experiment 
with  better  instrumentation  is  required  to  obtain  data  closer  to  the  boundaries  to  be  able  to 
finally  discern  the  cause  of  the  non-zero  velocity  at  the  guidewire  wall  and  at  the  entrance 
of  the  stenosis. 

COMSOL  proved  itself  to  be  capable  for  both  a  time  dependent  study  of  the 
model  as  well  as  for  the  application  of  pulsatile  flow.  Under  some  circumstances  an 
adaptive  mesh,  which  significantly  increased  processing  time,  was  necessary  in  order  to 
obtain  accurate  results. 

In  order  to  continue  developing  the  model  for  cardiovascular  purposes  an  accurate 
range  of  heart  beat  models  will  need  to  be  developed  in  order  to  ensure  that  the  system  is 
applicable  to  modeling  a  wide  range  of  patients.  This  will  require  extensive  research  into 
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both  the  ranges  of  velocity  profiles  that  vary  between  patients  as  well  as  determining  the 
varieties  of  profiles  based  on  the  location  of  the  stent  within  the  body.  Additionally, 
work  will  be  required  to  study  non-Newtonian  fluids  in  COMSOL  in  order  to  more 
accurately  model  blood  flow. 

Once  a  heartbeat  has  been  sufficiently  modeled  a  fluid  structure  interaction  will 
need  to  be  developed  to  begin  modeling  the  desired  sensor.  COMSOL  has  this  capability 
although  it  has  not  currently  been  tested  for  the  models  used  in  this  thesis.  COMSOL 
also  has  the  capability  to  model  piezoelectric  materials  that  can  be  used  to  create  a  sensor 
that  generates  its  own  power. 
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